'''
Created on May 25, 2011

@author: oabalbin
'''
from collections import defaultdict
def read_mut_genomic_coord(ifile):
    '''
    Reads file genomic oordinates from the hmgd database.
    coordinates are from the hg18
    '''
    ifile = open(ifile)
    i=0
    acc_id, chr, strand, start, end = 0,1,2,3,4
    genomic_coords=defaultdict()
    for l in ifile:
        fields=l.strip('\n').split('!')
        HMGD_ID, CHRM, STRAND, START, END = \
        fields[acc_id],fields[chr],fields[strand],fields[start],fields[end]
        #genomic_coords[HMGD_ID]=[CHRM, STRAND, START, END]  
        bed=['chr'+CHRM, START, END, HMGD_ID, '0',STRAND,'&&'] 
        genomic_coords[HMGD_ID]=bed
        #print HMGD_ID
        #print ",".join(bed).replace(',','\t')

    return genomic_coords

def read_allmut_table(ifile, ofile, genomic_coords=None):
    '''
    '''
    ifile = open(ifile)
    ofile_bed=open(ofile+'.bed','w')
    ofile_pos=open(ofile+'.pos.base','w')
    ofile = open(ofile,'w')
    #i=0
    acc_id, mut_type  = 16,18
    single_mutation='M'
    regular_field_length=19
    #insertion, deletion, others
    hmgd_muts=defaultdict()
    
    
    for l in ifile:
        fields=l.strip('\n').split('!')
        if len(fields) < regular_field_length:
            print fields 
            continue
            
        if fields[mut_type] == single_mutation:
            if genomic_coords is not None:
                hmgd_id=fields[acc_id]
                try:
                    bed=genomic_coords[hmgd_id]
                    rec=bed+fields
                    hmgd_muts[hmgd_id]=rec
                    base_reference, base_alternative = fields[7].split('>')[0][-1],fields[7].split('>')[1]
                    ofile.write(",".join(rec).replace(',','\t')+'\n')
                    bed[2]=str(int(bed[2])+1)
                    ofile_bed.write(",".join(bed).replace(',','\t')+'\n')
                    #ofile_pos.write(bed[0]+'@'+bed[1]+'|%s>%s'%(base_reference,base_alternative)+',')
                    ofile_pos.write(bed[0]+'@'+bed[1]+',')
                except KeyError:
                    pass
                    
    return hmgd_muts

        
ifile = '/data/databases/HMGD/genomic_coords.txt'
mfile='/data/databases/HMGD/allmut.txt'
ofile='/data/databases/HMGD/hmgd_single_mutations.txt'
genomic_coords=read_mut_genomic_coord(ifile)
read_allmut_table(mfile, ofile, genomic_coords)
